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^ ! Abstract 



The Markov chain Monte Carlo methods offer practical procedures for detecting signals char- 
acterized by a large number of parameters and under conditions of low signal-to-noise ratio. We 
present a Metropolis-Hastings algorithm capable of inferring the spin and orientation parameters of 
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I. INTRODUCTION 



The world-wide network of laser interferometric gravitational wave detectors has begun 
to ac qui re seientineal.y signified data fl Q fl Q and rapid* bating neutron stars are 
an important potential source of signals (we will reserve the term 'pulsar' to refer to the 
observed pulsating radio sources). Although a spinning spherically symmetric neutron star 
will not produce gravitational waves, a number of mechanisms have been proposed that are 
capable of producing quasi-periodic gravitational waves from biaxial or triaxial neutron stars 
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Any gravitational waves from these neutron stars will likely be seen at Earth as weak 
continuous wave signals. 

The data analysis task of identifying such a signal in the output of a laser interferometer 
is challenging and difficult, both because of the weakness of the signal and because its time 
evolution is characterised by a relatively large number of parameters. Radio observations 
can provide the sky location, rotation frequency and spindown rate of known pulsars, but 
the problem of looking for unknown (or poorly parameterised) neutron star sources is signif- 
icantly more challenging. SN1987A is a good example of a poorly parameterised source for 
which the sky location in approximately known but also for which there is a large uncertainty 
in the frequency and spindown parameters of the putative neutron star 

Much work has already gone into all-sky hierarchical methods for searching for continu- 
ous gravitational waves 8|, |9|. Here we address the specific problem of a 'fuzzy' parameter 
space search, in which a restricted volume of the space needs to be thoroughly investigated. 
We take a Bayesian approach to this problem and use Markov chain Monte Carlo (MCMC) 
techniques which have been shown to be especially suited to similar problems invol ving nu- 
merous parameters [10]. In particular, the Metropolis-Hastings (MH) algorithm [ljj, [l2| has 
jeen used for estimating cosmological parameters from cosmic microwave background data 
14 , ^| , and the applicability of the MH routine has been demonstrated in estimating 
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astrophysical parameters for gravitational wave signals from coalescing compact binary sys- 




17]. MCMC methods have also provided Bayesian inference for noisy and chaotic 
19]. 



Here we demonstrate that a MH algorithm also offers great promise for estimating neu- 
tron star parameters from their continuous gravitational wave signals. This works builds 
on the development (by two of us) of an end-to-end robust Bayesian method of searching 
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for periodic signals in gravitational wave interferometer data [2Cj, summarized in Sec. |TTJ 
Starting with this Bayesian approach we apply a similar MH routine to that used in 
The description of the Bayesian MH method is described in Sec. 11111 In Sec. II VI we present 
the results of this study, using synthesized data, for four and five parameter problems. We 
believe that this method offers great hope for signal extraction as more parameters are 
included, and this point is discussed in Sec. El 

II. SIGNAL CHARACTERISTICS 

We will initially consider searching for signals from known radio pulsars, and then expand 
the method to account for an uncertainty in the frequency of the gravitational wave signal. 
As gravitational waves from pulsars are certainly weak at Earth, long integration periods 
are required to extract the signal, and we must take account of the antenna patterns of the 
detectors and the Doppler shift due to the motion of the Earth. 



As in the previous study 
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21( we consider the signal expected from a non-precessing 
triaxial neutron star. The gravitational wave signal from such an object is at twice its 
rotation frequency, / s = 2/ r , and we characterise the amplitudes of each polarization with 
overall strain factor, ho. The measured gravitational wave signal will also depend on the 
polarisation antenna patterns of the detector F x>+ giving a signal 

s(t) = -F + (t;ip)h (l + cos 2 i)cos^(t) + F x (t; ip)h Q cos l sin #(£), (1) 

where ip is the polarization angle of the radiation (which depends on the position angle of 
the spin axis in the plane of the sky) and i is the inclination of the pulsar with respect to 
the line-of-sight. 

Using a simple slowdown model, the phase evolution of the signal can be usefully param- 
eterised as 



= 0o + 2tt 

where 



f s (T - T ) + i/„(T - T ) 2 + ~/ s (T - T ) 3 



(2) 



T = t + 5t = t + ^ + A T . (3) 



Here, T is the time of arrival of the signal at the solar system barycenter, 0o is the phase 
of the signal at a fiducial time T , r is the position of the detector with regard to the solar 



system barycenter, n is a unit vector in the direction of the pulsar, c is the speed of light, 
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and At contains the relativistic corrections to the arrival time 

The signal is heterodyned by multiplying the data by exp(— i^?(t)) so that the only time 
varying quantity remaining is the antenna pattern of the interferometer (which varies over 
the day). For convenience, the result is low-pass filtered and resampled. We are left with 
a simple model with four unknown parameters: the overall amplitude of the gravitational 
wave signal (ho), its polarization angle (tp), its phase at time T (0 O ) , and the angle between 
the spin axis of the pulsar and the line-of-sight (l). 

nn 

A detailed description of the heterodyning procedure is presented elsewhere (20, |21| ; here 
we just provide a summary of this standard technique. The raw signal, s(t), is centered near 
twice the rotation frequency of the pulsar, but is Doppler modulated due to the motion of 
the Earth and the orbit of the pulsar if it is in a binary system. The modulation bandwidth 
is typically 10 4 times less than the detector bandwidth, so one can greatly reduce the 
effective data rate by extracting this band and shifting it to zero frequency. In its standard 
form the result is one binned data point, Bk, every minute, containing all the relevant 
information from the original time series but at only 2 x 10~ 6 the original data rate. If the 
phase evolution has been correctly accounted for at this heterodyning stage then the only 
time-varying component left in the signal will be the effect of the antenna pattern of the 
interferometer, as its geometry with respect to the neutron star varies with Earth rotation. 
Any small error, A/, in the heterodyne frequency will cause the signal to oscillate at A/, 
and for the second part of our study we have Af as our fifth parameter. For both these 
studies we estimate the noise variance, a\, in the bin values, Bk, from the sample variance 
of the contributing data. It is assumed that the noise is stationary over the 60s of data 
contributing to each bin. 



III. THE METROPOLIS-HASTINGS ALGORITHM 



This section presents a brief review of the Bayesian MH approach to parameter estima- 
tion. Comprehensive descriptions of MCMC methods and the MH algorithm can be found 
elsewhere HHIlJ- 

We will denote the output from the above heterodyning procedure as {Bk}, with joint 
probability distribution function (pdf) denoted by p({Bk}\a) conditional on unobserved 
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parameters a = (aj., . . . , a^). The pdf p({B k }\a) is referred to as the likelihood and regarded 
as a function of the parameters a. The parameters of interest for our four parameter study 
are a = (h Q , ip, (fi , t), while for the five parameter study they are a = (h , i/j, i, Af). 
From Eq. the (now complex) heterodyned signal is 

y(t k ;a) = ^F + (t k ; t/;)h (l + cos 2 t)e^ - ~ F x (t k ;^)h cos te^ , (4) 

and the binning procedure should, by the central limit theorem, give the noise a near- 
gaussian probability density characterized by a variance a\ for the kth bin. The likelihood 
that the data in this bin, taken at time t k , is consistent with the above model is 



\B k -y(t k ;a)\ 



p{B k \a) oc exp ( ) , (5) 



and the joint likelihood that the data in all the bins (taken as independent) are consistent 
with a particular set of model parameters is 

p({B t }|a)ocnexp( " |Bt " 2 y ;a)|2 ). (6) 

Bayesian inference requires the specification of a prior pdf for a, p(a), that quantifies the 
researcher's pre-experimental knowledge about a. The phase and polarisation priors are flat 
in their space, and are set uniform for 0o over [0, ir], and for ip over [— 7r/4, vr/4]. The prior 
for l is uniform in cost over [—1, 1], corresponding to a uniform prior per unit solid angle of 
pulsar orientation. Finally, in the present study we take a prior for ho that is uniform for 
< h < 1000 (in our normalized units for which a k — 1), and zero for all other values. 

Using Bayes' theorem, the post-experimental knowledge of a is expressed by the posterior 
pdf of a: 

iKoHa.H- ^gW'-^ jK.MWI.), P) 

where p({B k }) = J p({B k }\a)p(a) da is the marginal pdf of {B k } which can be regarded as 
a normalizing constant as it is independent of a. The posterior pdf is thus proportional to 
the product of prior and likelihood. 

The marginal posterior distribution for parameter is the integral of the joint posterior 
pdf over all other components of a other than a^, i.e., 

p{ai\{B k }) = J ... J p(a\{B k }) dai . . . da^i da m . . . da d , (8) 
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and contains all the analysis has to say about the value of dj alone. However it is often 
useful to summarise this in a single 'point estimate' of dj using, for example, the posterior 
mean: 

(a*) = J aip(ai\{B k }) ddi. (9) 

Calculating the normalization constant p({Bk}) and calculating each marginal posterior pdf 
requires difficult d- and d — 1 dimensional integrations, respectively, that can be evaluated 
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17|. Rather than sampling di- 



usi „ g a sampling _ h and MCMC m e thods U 
rectly from p(a\{Bk}), a sample from a Markov chain is generated which has p(a\{Bk}) as 
its equilibrium distribution. Thus, after running the Markov chain for a certain 'burn-in' 
period, these (correlated) samples can be regarded as samples from the limiting distribu- 
tion, provided that the Markov chain has reached convergence. Despite their correlations, 
the ergodic theorem guarantees that the sample average is still a consistent estimate of the 
posterior mean Eq. Q 23]. 



11 



121. The MH 



The specific MCMC technique used for this study was the MH algorithm 
algorithm generates a sample from the target pdf p(a\{Bk}) using a technique that is similar 
to the well-known simulation technique of rejection sampling. A candidate is generated 
from an auxiliary pdf and then accepted or rejected with some probability. Starting with 
an arbitrary initial state clq, at time n a new candidate a' is generated from the candidate 
generating pdf, q(a\a n ), which can depend on the current state a n of the Markov chain. This 
new candidate a 1 is accepted with a certain acceptance probability a(a'\a n ), also depending 
on the current state a n , given by 

//in • / p{a')p({B k }\a')q(a n \a') \ 
ala \a n ) = mm < — — — - - — — — — — —, 1 > . 10) 

\p(a n )p({B k }\a n )q(a'\a n y J 1 ; 

For good efficiency a multivariate normal distribution centered at the current state a n is 
used for q(a'\a n ). This then implies that if the posterior probability at a' is larger than at 
the current state a n , the proposed step to a' is always accepted. However, if the step is 
in a direction of lower posterior probability, then this step is accepted only with a certain 
probability given by the ratio of the posterior pdfs (since our multivariate normal generating 
function is symmetric in a' and a n and therefore cancels out). If the candidate is accepted, 
the next state of the Markov chain is a n+ i = a', otherwise the chain does not move, i.e. 

The steps of the MH algorithm are therefore: 



Step 0: Start with an arbitrary value do 
Step n+ 1: Generate a' from g(a|a n ) and u from £7(0, 1) 
If u < a(a'\a n ) set a n+1 = a' (acceptance) 
If u > a(a'\a n ) set a n+ i = a n (rejection). 
The efficiency of the MH algorithm depends heavily on the choice of the proposal density. 
The closer the proposal is to the target distribution, the faster convergence will be accom- 
plished. This link between the closeness of the proposal to stationary distribution and speed 
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of convergence has also been substantiated by Holden 
dynamically altered the proposal distribution based on information from the chain's history. 
The approach, called pilot adaptation, is to perform a separate pilot run to gain insight 
about the target density and then tune the proposal accordingly for the successive runs. 
Such adaptation can be iterated but allowing it infinitely often will destroy the Markovian 
property of the chain and thereby often compromise the stationarity of the chain and the 
consistency of sample path averages ( 2^; see for an example). 

Based on the central limit theorem, the posterior pdf should be well approximated by 
a multivariate normal distribution with mean equal to the posterior mode and covariance 
matrix equal to minus the Hessian evaluated at the posterior mode. Thus, we use a mul- 
tivariate normal distribution for the proposal density q(a\a n ). As the mode is unknown, 
we try to make use of pilot samples to estimate its covariance matrix. When we initially 
run the MH algorithm, we sample candidate parameters from a normal distribution with 
covariance matrix equal to the identity matrix and centered around the current state. After 
the completion of this pilot run we use the empirical covariance matrix of the sample as 
covariance matrix of the multivariate normal proposal density, again with mean equal to the 
current state. 



IV. RESULTS 

In the first part of our study we reproduced the results presented in |2j| where the four 
unknown parameters were ho, i, ip, and 0o- The signal s{t) was synthesized assuming a source 
at RA = 4 h 41 m 54 s and dec = 18° 23' 32", as would be seen by the LIGO-Livingston interfer- 
ometer. This was then added to white gaussian noise, n(t), which is a good approximation 
to the detector noise in our band. Our normalized data had a noise variance of o~\ — 1 for 
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each sample, and the amplitude of the signal used in our test runs was varied in the range 
ho = 0.0 to 10.0. We were able to detect signals for h > 0.1. The length of the data set 
corresponded to 14400 samples or 10 days of data at a rate of one sample per minute (which 
was the rate used for the LIGO/GEO SI analysis described in 2l|). Although we will work 
with strains normalised to = 1, the results can be cast into a more conventional form by 
multiplying cjj and ho by (S , / t /60) 1//2 , where (Sh) 1 ^ 2 is the strain noise spectral density of the 
detector at the frequency of interest, in Hz -1 . 

An example of the MH routine output is shown in Fig. ^ Displayed are the trace plots 
and the kernel densities (posterior pdfs). For this example the program ran for 10 6 iterations. 
The first 10 5 iterations were discarded as the burn-in. Short-term correlations in the chain 
were eliminated by 'thinning' the remaining terms; we kept every 250 th item in the chain. 
The true parameter values for this run were ho = 5.0, ip = 0.4, <po = 1.0 and i = 0.5 
(cos i = 0.88). In the example displayed in Fig. [T]the MCMC yielded mean values and 95% 
posterior probability intervals of h = 4.9 (4.43 to 5.50), if; = 0.02 (-0.68 to 0.69), O = 1-34 
(0.71 to 2.08), and cost = 0.90 (0.79 to 0.99). The 95% posterior probability interval is 
specified by the 2.5% and 97.5% percentile of p(ai\{Bk}). In Fig.Elwe display the estimated 
posterior pdf of h on an expanded scale, along with the real and estimated value for ho- 
lt is crucial that our algorithm is sensitive to the true value of the gravitational wave 
amplitude, h , even under conditions of relatively low signal-to-noise ratio, and Fig. El shows 
injected ho values versus their values inferred by the MH routine. The error bars correspond 
to the 95% posterior probability interval, i.e. the lower and upper bound are specified by 
the 2.5% and 97.5% percentile of p(a>i\{Bk}). The algorithm clearly is successful in finding 
and estimating h . While the error bars increase as the signal gets larger, the relative error 
Nio/ho does diminish as ho increases. The fact that the 95% posterior probability interval 
grows with ho for constant noise level would seem to be counterintuitive. In addition, the 
widths of the posterior probability distributions for ho are larger than would be naively 
expected from a search for a simple periodic signal. The reason is that these error bars 
represent the uncertainty in the parameter rather than just the level of the noise, and this is 
affected both by the noise level and the posterior covariance between all of the parameters. 
The MCMC technique also allows one to calculate cross-correlation coefficients from the 
Markov chains of the parameters, and the value between h and cos i in all of our runs was 
~ — 0.95. As a result the data are consistent with a relatively broad range of combinations of 
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FIG. 1: Trace plot (left) and MCMC-estimated posterior pdfs (right) for the pulsar parameters 
ho, if}, 4>o and cosi. In this example the true parameters were ho = 5.0, (fio = 1.0, ip = 0.4, and 
i = 0.5 implying cos i = 0.88. 

the two parameters, making their individual values rather uncertain here - an effect evident 
from Eq. 

The effect of the other unknown parameters (particularly i) on the posterior pdf for ho 
can be clearly shown by repeating the analysis for Fig. |3] but with ho as the only unknown, 
namely, all of the other parameters set to their actual values in the MCMC routine. Under 
these circumstances the widths of all 95% posterior probability intervals are 0.116, indepen- 
dent of the value of h . Comprehensive analyses have investigated detection statistics for a 
periodic signal in a gravity wave detector [27[. However, these statistics are concerned only 
with the amplitude of the periodic signal, and not with parameter estimation (as described 
above). If we write Eq. ^as s(t) = Acos(^ + $) (with A being the periodic signal ampli- 
tude and $ a phase term) then the detection statistic of j^] would apply to finding a signal 
amplitude A in the presence of the detector noise. In terms of Eq. ^ the amplitude of the 
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FIG. 2: An expanded view of the estimated posterior pdf based on the MCMC sample for param- 
eter ho- The vertical solid line shows the posterior mean of ho = 4.9, while the vertical dotted line 
marks the true parameter value of ho = 5.0. 
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2 4 6 8 10 

true ho 

FIG. 3: The posterior mean based on the MCMC sample for the gravitational wave amplitude 
parameter ho versus the actual value of ho used in synthesizing the data. The error bars correspond 
to lower and upper bounds at the 2.5% and 97.5% percentiles of the posterior pdf. The solid line 
has a slope of 1. The calculations were performed over 14 000 data points, each with noise variance 
of a\ = 1. 

periodic signal would be 

A = {[F + (t; ^)A (1 + cos 2 i)/2} 2 + [F x (t; ^)h cos if j 1 / 2 . (11) 

It is clear that A has a complicated dependence on ho and cos i. We will never know, a priori, 
the value of all the pulsar parameters. Our study here is about parameter estimation, and 
not knowing the values of all the pulsar parameters ultimately increases the width in the 
posterior pdf for the gravity wave magnitude ho- 

As the magnitudes of the signals are diminished there comes a point when one is no longer 
able to confidently claim a detection. This threshold is somewhat arbitrary, and dependent 
on the statistics and interpretation. In the study presented here we claim that a signal is 
detected when the ho = point is more that two standard deviations from the mean value of 
the MCMC generated posterior pdf for ho- For the synthesized signals we investigated this 
corresponded to a threshold for detection of h = 0.1; in this case the measured mean of the 
posterior pdf for ho was 2.1 standard deviations away from zero. For an initial detection of 
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gravitational radiation it is likely that the scientific community will demand a significantly 
larger signal-to-noise ratio. However, the performance of the MCMC routine is still very 
good for these relatively low signal levels. 

Although 10 6 Monte Carlo iterations were used in this study (taking 1 d on a 1 GHz 
processor) adequate distributions can be generated from 10 5 iterations after the burn-in, so 
good results can be achieved after just a few hours. In fact the marginalisations discussed 
above can be tackled more quickly using simple summing methods as performed by j^, 
and the result of a comparison of the two is shown in Fig. |3J The great advantages of the 
MCMC method for us is its demonstrated ability to deal with problems that have a large 

, where other numerical integration techniques (such as employed 
by [si]]) are not feasible. The ultimate goal of our research is to expand this pulsar parameter 
estimation work to include more parameters. The next step in increasing the complexity 
of the pulsar signal is to consider potential sources of known location, but with unknown 
rotation frequency. In order to start this investigation we added a new parameter, the 
uncertainty in the frequency of the source, Af. In this example the exact value of the 
pulsar's gravitational wave signal is uncertain to within 1/60 Hz. In the study we present 
here there is a difference, Af, between the gravitational wave signal frequency and the 
heterodyne frequency. The addition of this new parameter did not significantly increase the 
rate at which the code ran, but did (by about 20%) increase the length of the burn-in time. 
If one wanted to increase this frequency range to 5 Hz then this could be done by running 
the MCMC code on 300 processors, with each run differing in center frequency by 1/60 Hz. 
The Markov chain using the correct frequency would converge, while the other 299 chains 
would not. This will be a future research project for us. 

In our MH code we used a uniform prior for the uncertainty in the frequency, Af, over 
±0.016 67 Hz. The injection parameters used were ip = 0-4, 4> = 1-0, Af = 0.007812 5, 
and i = 0.5 (cos i = 0.88). ho was again injected with a number of values between 0.25 
and 10.0. In Fig. 03 we show sample trace plots and posterior pdfs for Af and ho when the 
injected value of ho was 1.0. For this example the MCMC algorithm yielded mean values 
and 95% posterior probability intervals of h = 1.02 (0.86 to 1.26) and Af = 0.007812 497 
(0.007812 480 to 0.007812 515). The frequency pdf is quite narrow, which was responsible 
for the increase in the burn-in time as the Markov chain must find this narrow region of 
parameter space. In Fig. |H1 we display the estimate for the gravitational wave amplitude 
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FIG. 4: The posterior mean based on the MCMC sample for the gravitational wave amp litude 
parameter ho (dotted line), along with the that produced via the method presented in [20j (solid 
line) In this example the true value was ho = 0.5, while the other true parameter values were 
$ = 0.4, <j) = 1.0, and t = 0.5. 

(ho) predicted by the five parameter MH routine versus the actual ho. In Fig. [7] we display 
the estimate for the difference in frequency Af predicted by the five parameter MH routine 
versus the injected ho- 



V. DISCUSSION 



Recent applications of MCMC techniques have provided a Bayesian approach to estimat- 
ing parameters in a number of physical situations. These include cosmological parameter 
estimation from cosmic microwave background data 



14, 



151 ] - estimating astro phy sical 

parameters for gravitational wave signals from coalescing compact binary systems 16L Il7|. 
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and parameter estimation of a chaotic system in the presence of noise |18l. Il9j. An all sky 
survey for periodic gravitational waves from neutron stars must explore a very large param- 
eter, and this has partially been addressed in Q|. Generically, the signal from a neutron 
star in a binary system will be characterized by at least 13 parameters. It is our hope that 
MCMC techniques will prove fruitful in dealing with these complex signals. 
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FIG. 5: Trace plot (left) and posterior pdfs (right) for the pulsar parameters ho and Af. In 
this example from the five parameter problem the true values for these critical parameters were 
h = 1.0 and Af = 0.0078125. 

In this paper we have demonstrated that the success of MH routine for the five parameter 
problem: ho, if), <f>o, i and Af. Our longer term plans are to account for other parameters, 
such as spindown rate, pulsar wobble, and possibly location of the signal in the sky. This 
research is currently in progress. 
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FIG. 6: The posterior mean based on the MCMC sample for the gravitational wave amplitude 
parameter ho versus the actual value of ho used in synthesizing the data. This example is from the 
five parameter problem. The error bars correspond to the lower and upper bound being specified 
by the 2.5% and 97.5% percentiles of the posterior pdf. The solid line has a slope of 1. 
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